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Analytical study of electrostatic ion beam traps 
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The use of electrostatic ion beam traps require to set many potentials on the electrodes (ten in 
our case), making the tuning much more difficult than with quadrupole traps. In order to obtain 
the best trapping conditions, an analytical formula giving the electrostatic potential inside the trap 
is required. In this paper, we present a general method to calculate the analytical expression of the 
electrostatic potential in any axisymmetric set of electrodes. We use conformal mapping to simplify 
the geometry of the boundary. The calculation is then performed in a space of simple geometry. 
We show that this method, providing good accuracy, allows to obtain the potential on the axis as 
an analytic function of the potentials applied to the electrodes, thus leading to fast, accurate and 
efficient calculations. We conclude by presenting stability maps depending on the potentials that 
enabled us to find the good trapping conditions for O^^ at much higher energies than what has 
been achieved until now. 

PACS numbers: Valid PACS appear here 



I. INTRODUCTION 

In the last few years a variety of electrostatic devices 
for storing and handling low energy ion beams have been 
designed and operated. Electrostatic storage rings [l|-l3|, 
cone traps ^.^rhitraps [5], and electrostatic ion beam 
traps (EIBT)[6|, [3 are now used to study atomic and 
molecular metastable states or molecular fragmentation, 
photodissociation or mass spectrometry (see, e.g., Ref. 
Q for a review). The design and study of these instru- 
ments relies nowadays mainly on computer simulations. 
An EIBT, as designed by D. Zajfman and collaborators 
[6',"9'], and independently by W.H. Benner [7], is a purely 
electrostatic trap composed of two electrostatic mirror- 
Einzel lens combinations, as represented on Fig. [TJ This 
trap has many interesting features [iO| : on one hand, it 
offers trapping of energetic particles (keV) in a well defi- 
ned direction and on the other hand it is small, relatively 
inexpensive and has a field-free region where ions move 
freely and where measurements can easily be performed. 
It can also be used as a moderate-resolution mass spec- 
trometer [ll|. 

In this paper, our aim is to provide a method enabling 
the determination of analytical solution to the electro- 
static potential in any axially symmetric configuration 
using the elegant method of conformal mapping. We will 
present the method on the EIBT, but it can be used on 
other sets of electrodes. 

The main reason for quadrupole traps' extensive use 
in precision experiments, lies in the fact that their fields 
can easily be described by an analytical formula. It en- 
ables a deeper understanding of many subtle phenomena 
like frequency shifts due to space charge. The inventors of 



the ion trap resonator used either a matrix approach [ll| 
or numerical simulations [12|. The former are useful for 
acquiring a qualitative understanding, but cannot pro- 
vide detailed insight in the operating conditions, while 
the later may not provide enough numerical accuracy to 
follow the particles during several tens of thousands of 
oscillations inside the trap, or be too time and resources 
consuming to explore many different potential configu- 
rations. The tracking of particles in the EIBT is rather 
difficult because of the combination of a long free-flight 
zone between the mirrors and of two areas in which the 
particles slow down, stop and reverse their course on very 
short distances, while being subjected to strong, rapidly 
varying, electrostatic fields. The simulation time becomes 
rapidly a limitation when many potential configurations 
must be studied while looking for new operating condi- 
tions. Instead of two tuning parameters as in a Paul trap, 
we have to fix the potentials of five electrodes in a sym- 
metric configuration (ten when each side of the trap is set 
differently) . The space to explore is therefore too large for 
numerical simulations in which the field is determined by 
usual finite element methods, where one has to make a 
different calculation for each set of parameters. Moreover 
numerical errors accumulate and perturb the trajectory 
of the particle over long trapping times. In the sequel, we 
will show how to find a formula, depending only on these 
parameters, which is able to give the electrostatic poten- 
tial with good accuracy. We will also show some practical 
applications. 

Before we present our method, we will just review two 
ways of calculating electrostatic potentials and explain 
why they are difficult to use in our case. 

Green's functions often yield analytical results because 
they allow reduction of the solution of the Dirichlet pro- 
blem to the calculation of the following integral : 
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Figure 1: (Color online) Overview of the EIBT. Five poten- 
tials are applied to the electrodes, the other are grounded. 
The injection of the bunch of ions is performed when all the 
electrodes on one side of the trap are grounded. The potentials 
are raised before the bunch has time to come back. This trap 
can be used to make met ast able state lifetime measurements, 
hence the photomultiplier tube. 



where Mq is the point where the potential is evaluated, eo 
is the vacuum permeability, n parametrizes the direction 
orthogonal to the surface, Us{M) is a given potential dis- 
tribution over the surface S and G is a Green's function. 
However, Green's functions are only known for simple 
geometries, which limits the analytical approach. 

A very interesting method called quasi- Greens' func- 
tions has been developed in [i3|. This method is based 
on the division of a complicated geometry into different 
simple shapes. The main drawback is that the final ex- 
pression is given as an infinite sum whose coefficients 
have no closed form. In practice, the given expression, 
although analytical, is much more complicated than the 
one we will present in this article. 

Another technique, is the charge ring method |14J17|. 
The integral form of the Poisson equation is applied to 
N rings representing the geometry. If N is large enough, 
we get a set of linear equations ^ = AQ , where <l> is the 
potentials applied to the rings, Q is the charge induced 
on each ring and A is a matrix depending only on the 
geometry. Once the inverse of A has been determined, 
the charge of each ring is know and the potential at a 
point r that is not on the boundary is given by : 
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where Si is the area of ring i. Hundreds of rings a usually 
provide a accuracy of the on-axis potential of the order 
of 10~^ (one order of magnitude better than the method 
proposed here, see Sec. lIII A|) . However, it is at least one 
order of magnitude slower to compute : even though the 
inverse of A is calculated only once for a given geome- 
try, one has to evaluate numerically N integrals, which 
is much longer than to evaluate common mathematical 
functions. We have implemented both methods on the 
same geometry and the conformal method was 64 times 
faster than the ring method. The choice between those 



two techniques will rely on the need to improve accuracy 
or speed. 

This article is organized as follows : in Sec. |TT1 we 
present an approximate method to obtain the analytical 
potential of an axisymmetric set of electrodes, having the 
same radius. In Sec. IIIII we explain how to use Schwarz- 
Christoffel method to alter the metric and obtain a set 
of electrodes with the same radius and in Sec. IIII Al we 
solve the problem in this new space using the method of 
section [TTl Section IIIIBI contains a summary of the key 
steps of the whole method as well as a discussion on the 
improvement of the accuracy. Finally, in Section IIVI we 
show that the dynamics of the ions in the EIBT is go- 
verned by an Hill's equation and we present a stability 
map showing what experimental parameters lead to an 
efficient trapping. 



II. SEPARATION OF VARIABLES AND 
BERTRAM'S METHOD 

We start from the Laplace equation in cylindrical co- 
ordinates : 
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where V = F (r, z) is the potential at radius r from the 
axis and at a distance z from the center of the trap. 
Using the method of separation of variables, i.e., assu- 
ming V(r, z) = R{r)Z{z)^ we obtain [18] : 
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where /c is a real constant. ([2]) is a particular case of the 
general Bessel equation [19], whose solution is a Bessel 
function of first kind : Jo{kr). We can use the general 
solution given by a Fourier-Bessel series of the form : 
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enabling us to take into account the boundary conditions. 
Thus, if the potential at some radius R is known as a 
function of z^ then a{k) can be found by means of the 
Fourier transform : 
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Given that any solution of the Laplace equation in a 
cylindrical symmetry is also of the general form : 
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it is sufficient to determine V(0,z) along the axis. 
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Figure 3: (Color online) A section of one side of the trap : the 
crosshatched areas are the electrodes, the continuous lines are 
the equipotentials, the domain A is surrounded by a dashed 
border and the circled points represent the Zi points, which 
define our boundary conditions. 



Figure 2: (Color online) The line represents the analytical 
function and the dots is the numerical solution achieved with 
the finite elements software COMSOL Multiphysics® (top). 
The difference between the two previous curves (bottom). 
In this case, we applied Bertram's method replacing R by 
a function of z. Even though the function is a polynomial 
going smoothly from Rz to R, we see that the error is large 
(-258V) near z = —0.15 where the transition occurs. Here 
{Vi = 4513 y, V2 = 4836 V, V3 = 3112 V, Va = 1642 V, 14 = 
,3941 V}. 



The previous method is rather simple and efficient to 
find the potential. Its only limitation is the need to have 
identical radii for all electrodes. In the next section, we 
will show how to use conformal mapping to place our- 
selves in the space where the borders have a constant 
radius. 



III. CONFORMAL MAPPING 



Following Bertram [20|, we assume that if we know 
the potential V{R^ Q at a distance R from the axis then, 
the potential on the axis is well approximated by the 
formula : 
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where the constant uj = 4Aq = 1.3152, and Aq is the first 
coefficient of the following development [20| : 
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In our case, if all the electrodes had the same ra- 
dius^ we could use directly this method, taking V{R, Q 
as a piecewise linear function of the set of potentials 
{"^1, ^2, "^3, V4, V2}. However, since the radius of the first 
four electrodes is R = 8mm, different from the Einzel 
electrodes where Rz = 13mm, this method does not work 
in the area where the radius changes. We tried to in- 
troduce a smooth function R{z) in the integral, but the 
difference between those results and a finite element so- 
lution always shows a large discrepancy as illustrated on 
Fig.H 



Conformal mapping is widely used in applied physics 
and chemistry. One can cite the design of airfoils : the 
Joukowsky transformation [21] reduces the study of the 
laminar flow on a complicated profile to the much ea- 
sier study of a cyHnder in the transformed flow, or the 
study of diffusive flow at micro-ring electrodes in ana- 
lytical chemistry [22|, l23|. Applications in electrostatic 
potential determination are also known (see, e.g., [2J- 
[27]). However, all these works use conformal mapping in 
a geometry where one dimension can be considered in- 
finite. A section perpendicular to this infinite dimension 
is then mapped. In this paper, we show that conformal 
mapping can be used to find the electrostatic potential in 
a space limited axisymmetric geometry by mapping the 
plane parallel to the axis and rotating it. 

In order to use the results of [Til the first step is to find 
the holomorphic function that maps section A, described 
on Fig. [31 on a rectangle. Holomorphic functions are of 
great interest because angles are conserved under those 
transformations : equipotential lines stay orthogonal to 
field lines. 

We first use the Schwarz-Christoffel transformation 
[23] to map domain A, parametrized by z = x -\- iy onto 
the upper half plane where t = r -h ^5, as represented on 
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Figure 4: (Color online) The three working planes. The Z- 
plane is the physical plane corresponding to domain A on 
Fig. O The T-plane is an intermediate step due to the fact 
that the Schwarz-ChristofFel method always maps a polygon 
on the upper complex plane where Im{t) — s > 0. The W- 
plane is the calculation plane where we can apply the method 
of section [Til since it can be seen as the section of a constant 
radius cylinder. The arrows on the border represent the three 
mapping used to transform one plane to the other. 



Fig. 1 : 
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f{t)=Krt-\t-ay/\t-b)-'/^ 



so that : 



where Ki and K2 are constants to be determined. With 
an appropriate choice of the origin (f(0)=0), we have 
i^2 = 0. For i^i, we use the boundary conditions : going 
from point A to point E in the Z-plane (see FigH]) cor- 
responds to a large semicircle of radius p -^ +(X) and 
from to TT in the T-plane : 
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When taking the limit p ^ +00, this reduces to : 
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and so Ki = Rz/'t^- The second boundary condition, 
going from B to B' in the Z-plane, is expressed by inte- 
grating around BB' with p ^ and 6 going from to 
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Choosing a = 1 only fixes the origin in the T-plane, and 
it implies Vb = ^ . Finally, we make the substitution : 



f -1 



and the integration gives : 
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The mapping from the W-plane to the T-plane is much 
simpler : 

w = g{t) = Logit) 

and we can finally link the Z-plane to the W-plane by 
the following transformation : 
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= CM, 



The transformation z = C{w) gives a one-to-one mapping 
between points of the Z-plane and points of the W^-plane 
and it is a conformal transform as the composition of two 
conformal transforms. Since ( is not invertible, we have 
found a good approximating function defined on three 
intervals and inverse of which is : 
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where zq = 0.007mm and P{z) is a 6th-order polynomial 
connecting continuously the two asymptotic expansions. 
The points at ±zo represent the limit of validity of the 
two asymptotic expansions. 

There is one last thing to do before one can solve the 
problem in the W^-plane. Fig. [5] shows the image of the 
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placing in (j6]), we obtain : 



Figure 5: (Color online) This figure represents the image of 
the mapping C('^) applied to the rectangle defined as w — 
u -\- iv with —20 < 1^ < 8 and — tt < v < 0. The dots indicate 
the position of the Zi. 



rectangle defined diS w = u ^ iv with —20<u<8 and 
— TT < V < 0. On the vertical bold lines, u is constant 
and we will call these lines iso-u. We see that the iso- 
u are distorted near 0. Therefore, the distance between 
the points Wi = (~^{zi) is not the same as the distance 
between the points Zi : the metric is not conserved on the 
border. Consequently, the width of the electrodes near 
the point C in the VF-plane is not the physical width, 
which leads to errors in the potential. In order not to 
modify the metric on the border, one should first apply 
the function 



/3{z) = ac\z)-m) 
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to the points Zi so that their image have the same dis- 
tance between them in the W^-plane and in the Z-plane 
(physical plane). This function uses the fact that, on the 
axis, the metric is not modified from one plane to the 
other. It is the same idea as in [26| where a "space de- 
pendent diffusion coefficient" is used to account for the 
fact that real space is compressed/expanded unevenly to 
fit the VK-plane. In the sequel, we shall use the following 
notation : Zi = l3{zi). 



A. Solving in the W-plane 

Using C{w) we can now apply Bertram's method in 
the W-plane where the radius is constant. Following |2Q| . 
we assume that the potential varies linearly between the 
electrodes. We shall come back to this approximation in 
the last part of this section. We now use (|6]) in the W- 
plane where the radius is constant R = tt. On the border, 
we have : 
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with Ui = C~^{^i)i where Zi = /3{zi). Zi are the posi- 
tions of the points on Fig. [3l ]^ is the potential at Zi. 
Il{ui -^ Ui-^i) is a function equal to zero everywhere ex- 
cept between Ui and Ui-^i where it is equal to one. Re- 
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Now that we have Vw^ the potential along the axis in 
the Z-plane is given by : 

Vz{S)^z) = Vw{C\^)) 

as k; = C~^(^) varies from point B to point A. 

The error between the numerical solution and our re- 
sult is around 1%, which is enough for many applications. 
There are two main sources of errors. First, the calcula- 
tion does not take into account the field leaking at the 
front and at the rear of the set of electrodes. It can be seen 
on Fig. [3] that the potential is not exactly zero on the axis 
after the neutral electrodes preceding Vz and following 
Vi. This explains the two dark negative zones on Fig. [71 
Second, on the same figure, we see that the hypotheses 
we have made, concerning the variation of the potential 
between the electrodes, induce an error of about lOOV on 
the border. Yet, since this error is oscillating along each 
border, it compensates and the error on the axis is only 
about 20V (Fig. [6]). An attempt to enhance the potential 
at the borders is described in [29|. Finally, there is an 
error coming from the cylindrical term of Eq. ([!]) , which 
is not invariant under conformal mapping. However, it is 
always at least one order of magnitude smaller than the 
two previous error types. 

The best way to increase accuracy is to use a domain 
that takes into account the shape of the electrodes. Ins- 
tead of the mapping of Fig. [5l we could have used the 
one of Fig. [HI Besides showing that this method can be 
applied to comphcated geometries. Fig. [8] enables us to 
emphasize a crucial point : the mapping might not be 
analytic for a given geometry. However, since the geo- 
metry does not change, it is enough to calculate the in- 
verse map numerically one time, find the polynomial that 
fits this numerical solution on the axis and then, using 
Bertram's formula we have an analytic expression whose 
parameters are the electrodes' potentials. With this map- 
ping, the accuracy is approximately 0.1%. 




Figure 6: (Color online) The line represents the analytical 
function and the dots is the numerical solution (top). The 
difference between the two previous curves (bottom). We see 
that the error is around 1% and its smooth repartition shows 
that it does not arise from the varying radius. This difference 
is the consequence of the approximation we used to fix V{R, z) 
between the electrodes. Here {Vi = 4513 V, V2 = 4836 V, V2, = 
3112 V, Va = 1642 y, Vz =, 3941 V}. 



Solving numerically a Schwartz- Christoffel map will 
not be treated here because we chose a complete ana- 
lytical case to show all the details of the method. One 
could refer to [23, l30| for an extensive review on all the 
numerical techniques involved. Note that all the results 
presented in the sequel, use the simple case of Fig. [51 




Figure 7: The line represents the potential on the border used 
to calculate the potential on the axis and the dots is the nume- 
rical solution achieved with Comsol multiphysics (top). The 
difference between the two previous curves (bottom). Here 
{Vi = 4513 V, V2 = 4836 V, V3 = 3112 V, Va = 1642 V, 14 = 
,3941 V}. 
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Figure 8: This figure represents the image of a mapping ob- 
tained numerically. 



B. Summary of the method 



5. Use the formula ([TQ|) to obtain Vwiu^ —Pi). 



1. Write the the Schwartz- Christoffel transformation 
adapted to the geometry to obtain z = C('^)- For 
more details see (2^. 



6. The potential of one point z = x-\-Oxion the axis 
in the Z-plane is Vz{z^O) = Vw{C~^{^))' 



2. If the inverse w = ( ^{z) is not straightforward, it 
can be well approximated by a polynomial. 

3. Let Zi = Xi -\- iR{xi) be the points of the Z-plane 
defining the position of each side of the electrodes, 
one has to apply the function I3{z) = ({(~^{z) — i7r) 
to obtain the set of points Zi = C{C~^{^i) ~ ^^) 

4. Transpose the problem from the Z-plane to the 
W-plane where Wi = ({zi) = Ui -\- ivi. 



IV. APPLICATION : STABILITY MAP 

The main difficulty to tune an EIBT is due to the 
large number of parameters implied in its manipulation : 
five potential values on each side, the energy of the ions, 
their temperature and their charge-to-mass ratio. The 
designers of the trap used an optical model consisting 
of mirrors and lenses [10|, yet, since the focal length is 
not linked to the values of the potential, the behavior 
of the EIBT, given a set of potentials, is unpredictable. 
Until now, experimentalists had to use simulation soft- 
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Figure 9: (Color online) These curves show 1 — | A| as a func- 
tion of Vz (potential applied to the lens) in two distinct cases : 
(right curve, +) is the trapping efficiency of Ar'^ at 4.2keV 
with the same conditions as in i32]- The set of points is an 
reproduction of the experimental data in their Fig. 3 a. (left 
curve, o) represents the trapping efficiency of Ar~^ at 1.2keV 
with the same conditions as in 10]. The set of points is an 
reproduction of the experimental data in their Fig. 8. There 
is very good agreement between theory and experiment. 



ware like SIMION@ ^] to determine optimal trapping 
conditions [lO|. Since a finite element calculation has to 
be achieved each time a parameter is changed, this me- 
thod turns out to be fastidious. Worse, simulating the 
trajectory of ions going back and forth is much more dif- 
ficult than simulating a beam because the errors on the 
position accumulate and become comparable to the size 
of the trap. Even with the smallest step and a recursive 
method, SIMION® [31] was not able to enforce energy 
conservation after a few hundreds oscillations (the error 
was about lOOeV). 

From Eq. (J3j), we can limit ourselves to the second order. 
Higher order terms can be neglected as long as r < Smm^ 
which is the aperture of our trap. The potential in the 
trap is then given by 



<|)(r, z) = V{z) r 



1 ,£V{z) 



dz^ 



(11) 



where V{z) is the potential along the axis. It follows that 
the trajectory of one ion in the trap is determined by the 
following set of equations : 
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The second term of the right side of equation ([12]) is 
small compared to the first (for r < 8 mm) and can thus 
be neglected. We obtain the longitudinal motion of the 
ion z{t). Substituting z{t) in (fT3]) . we obtain an Hill's 
equation [33] : 

d^r f q d^V{z) 



r = 0, 



(14) 



the term in parentheses being a periodic function of per- 
iod T. This equation arose in the study of the moon's 
dynamic and the usual method to discuss the stability of 
its solution consists of calculating infinite determinants 
|34| . The principal matrix of this equation is : 



M{t) 






(15) 



where V^i(t,to) is the solution of (p!4|) with initial condi- 
tions ^pl{to^to) = 1 and '0i(to,to) = and ip2{t^to) with 
'^2(^0,^0) = and '^2(^0,^0) = 1- Liouville's formula j35j 
shows that : 



detM{t,to) 



(16) 



and therefore the characteristic equation of the mono- 
dromy matrix M{to -\- T) is given by [35| : 



2Ax + 1 = 



(17) 



where : 



A = Tr{M{to + T)) = 



iPi{to+T,to)+Mto+T,to) 



1-1 (^^^ 

Applying Floquet's theorem [36|, we know that if A^ > 1, 

one of the two solutions is unbound but if A^ < 1 there 

are two solutions : 
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= ±7* 



P±{t) 



(19) 



where p±{t + T) = p±{t)^ 7 = Im{^\og{x-^)) and 
x+ = A + VA^ — 1. In conclusion. Hill's equation is 
stable, and thus trapping can be observed, when | A| < 1. 
We now compare these theoretical results with expe- 
riment. 

Figure [9] shows a comparison with the results publi- 
shed by other groups using the same kind of trap. These 
curves show 1 — | A | as a function of Vz (potential applied 
to the lens). They show a perfect agreement with two 
independent groups ^32| , [10| . We also tried to reproduce 
the data of [33] using H2 at l.OkeV and a negative 
potential on the Einzel lens. We also predict three 
stability intervals, approximately at the same position, 
because we have not taken into account the differences 
in the geometry of their trap. However we notice that 
the method works also with negative potentials. Our 
method enables to plot these curves in less than a second 
whereas the cited authors had to make a SIMION(R) [31| 
simulation for each point. 

Figure [10] shows the stability map depending on two 
parameters Vi and Vz (the respective potentials of the 
rear electrode and of the Einzel lens). This is the equi- 
valent to the famous Ince-Strutt diagram used to tune 
quadrupolar traps [38|. The white dots indicate settings 
where trapping was experimentally observed. We fixed 
the value of Vi and scanned Vz . Trapped ions go through 
a ring at the center of the trap and induce a current. 
We then analyze this amplified current with a spectrum 
analyzer : if we see a peak corresponding to the oscilla- 
ting movement of the ions, we mark this position with a 
white dot. The radius of each dot is proportional to the 
trapping efficiency. 

There is a shift between theory and experiment, which 
was first thought to be caused by a technical problem : 
it is difficult to monitor with a good accuracy our power 
supplies up to 8kV, especially because they raise in a 
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Figure 10: (Color online) The contours show the iso-?7 values 
where 77 = |A| — 1. The stability region is defined by 77 < 0. 
The white dots indicate settings where trapping was experi- 
mentally observed. The radius of each dot is proportional to 
the trapping efficiency. The thick white line is the border of 
the three stability region : there are three stability zones mar- 
ked I, II and III corresponding to the three types of orbits of 
Fig-ttH We used O^^ at 5.2keV/charge and the potential set 
is : {VuV2 = 5850 V, V3 = 4150 V, V4 = 1650 V, V^}. 



and, in the region where the ions can move, the error 
is never greater than 1%. This analytical expression of 
the potential inside this trap provides users with a much 
more powerful tool to study and optimize this novel kind 
of trap. As an example, we study the stability of the trap 
and show that they agree with experiments gi ving a fast 
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few nanoseconds. However, a recent improvement of the 
model, taking space charge effects into account, seems to 
explain this shift. These results will be presented somew- 
here else. Nonetheless, we could not explain the absence 
of trapping on the upper part of the stability zone II. We 
tested to see if higher order terms of Eq. (J3]) could ac- 
count for this, without success. No trapping can be seen 
in zone III, but as shown in Fig. [TTJ this region corres- 
ponds to very peculiar closed trajectories, which seem 
more theoretical than observable. 



Figure 11: Trajectories of 0^+ at 5.2 kV (top) : Stability zone 
I : {^1 = 7500 V, V2 = 5850 V, V3 = 4150 V, V4 = 1650 V, V^ = 
3650 V} . (middle) : Stability zone II : {Vi = 6500 V, V2 - 
5850 V,V3 = 4150 V.V^ = 1650 V,Vz = 4950 V}. (bot- 
tom) : Stability zone III : {V^i = 5100 V, V2 = 5850 V, V3 = 
4150 V, V4 = 1650 y, Vz = 5400 V"}. 



and easy way to predict trapping parameters. 



V. CONCLUSION 

In this article, we have given a method to calculate ana- 
lytical solution to the Laplace equation in axially sym- 
metric devices. This method is very general and can be 
used for various parts constituting a beam Hue. We suc- 
cessfully applied it to the EIBT and showed that the 
electrostatic potential is given by a formula depending 
on the five electrodes potentials. The formula has been 
compared with many finite element calculation where the 
potentials have been changed on the whole achievable 
range (from OV to SOOOF for each of the five potentials) 
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